/* Copyright 2007-2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Patrick van der Velde
 *
 * This file is part of dnAnalytics.  dnAnalytics is licensed under the 
 * Microsoft Public License. See License.txt for a complete copy of the
 * license.
 */

using System;
using System.Collections.Generic;

namespace dnAnalytics.LinearAlgebra.Solvers.Preconditioners
{
    /// <summary>
    /// This class performs an Incomplete LU factorization with drop tolerance
    /// and partial pivoting. The drop tolerance indicates which additional entries
    /// will be dropped from the factorized LU matrices.
    /// </summary>
    /// <remarks>
    /// The ILUTP-Mem algorithm was taken from: <br/>
    /// ILUTP_Mem: a Space-Efficient Incomplete LU Preconditioner
    /// <br/>
    /// Tzu-Yi Chen, Department of Mathematics and Computer Science, <br/>
    /// Pomona College, Claremont CA 91711, USA <br/>
    /// Published in: <br/>
    /// Lecture Notes in Computer Science <br/>
    /// Volume 3046 / 2004 <br/>
    /// pp. 20 - 28 <br/>
    /// Algorithm is described in Section 2, page 22
    /// </remarks>
    /// <include file='../../../../examples/LinearAlgebra/Solvers/Preconditioners/Ilutp.xml' path='example'/> 
    public sealed class Ilutp : IPreConditioner
    {
        /// <summary>
        /// The default fill level.
        /// </summary>
        private const double sm_DefaultFillLevel = 200.0;

        /// <summary>
        /// The default drop tolerance.
        /// </summary>
        private const double sm_DefaultDropTolerance = 0.0001;

        /// <summary>
        /// Returns the default fill level.
        /// </summary>
        public static double DefaultFillLevel
        {
            get
            {
                return sm_DefaultFillLevel;
            }
        }

        /// <summary>
        /// Returns the default drop tolerance.
        /// </summary>
        public static double DefaultDropTolerance
        {
            get
            {
                return sm_DefaultDropTolerance;
            }
        }
                
        /// <summary>
        /// The decomposed upper triangular matrix.
        /// </summary>
        private SparseMatrix m_Upper;
        /// <summary>
        /// The decomposed lower triangular matrix.
        /// </summary>
        private SparseMatrix m_Lower;
        
        /// <summary>
        /// The array containing the pivot values.
        /// </summary>
        private int[] m_Pivots;
        /// <summary>
        /// The fill level.
        /// </summary>
        private double m_FillLevel = sm_DefaultFillLevel;
        /// <summary>
        /// The drop tolerance.
        /// </summary>
        private double m_DropTolerance = sm_DefaultDropTolerance;
        /// <summary>
        /// The pivot tolerance.
        /// </summary>
        private double m_PivotTolerance;

        /// <summary>
        /// Creates an instance of the <c>Ilutp</c> class with the default settings.
        /// </summary>
        public Ilutp()
        { }

        /// <summary>
        /// Creates an instance of the <c>Ilutp</c> class with the specified settings.
        /// </summary>
        /// <param name="fillLevel">
        /// The amount of fill that is allowed in the matrix. The value is a fraction of 
        /// the number of non-zero entries in the original matrix. Values should be positive.
        /// </param>
        /// <param name="dropTolerance">
        /// The absolute drop tolerance which indicates below what absolute value an entry 
        /// will be dropped from the matrix. A drop tolerance of 0.0 means that no values
        /// will be dropped. Values should always be positive.
        /// </param>
        /// <param name="pivotTolerance">
        /// The pivot tolerance which indicates at what level pivoting will take place. A
        /// value of 0.0 means that no pivoting will take place.
        /// </param>
        public Ilutp(double fillLevel, double dropTolerance, double pivotTolerance)
        {
            if (fillLevel < 0)
            {
                throw new ArgumentOutOfRangeException("fillLevel");
            }

            if (dropTolerance < 0)
            {
                throw new ArgumentOutOfRangeException("dropTolerance");
            }

            if (pivotTolerance < 0)
            {
                throw new ArgumentOutOfRangeException("pivotTolerance");
            }

            m_FillLevel = fillLevel;
            m_DropTolerance = dropTolerance;
            m_PivotTolerance = pivotTolerance;
        }
                        
        /// <summary>
        /// Gets or sets the amount of fill that is allowed in the matrix. The
        /// value is a fraction of the number of non-zero entries in the original
        /// matrix. The standard value is 200.
        /// </summary>
        /// <remarks>
        /// <para>
        /// Values should always be positive and can be higher than 1.0. A value lower
        /// than 1.0 means that the eventual preconditioner matrix will have fewer
        /// non-zero entries as the original matrix. A value higher than 1.0 means that
        /// the eventual preconditioner can have more non-zero values than the original 
        /// matrix.
        /// </para>
        /// <para>
        /// Note that any changes to the <b>FillLevel</b> after creating the preconditioner
        /// will invalidate the created preconditioner and will require a re-initialization of
        /// the preconditioner.
        /// </para>
        /// </remarks>
        /// <exception cref="ArgumentOutOfRangeException">Thrown if a negative value is provided.</exception>
        public double FillLevel
        {
            get
            {
                return m_FillLevel;
            }
            set
            {
                if (value < 0)
                {
                    throw new ArgumentOutOfRangeException("Value");
                }

                m_FillLevel = value;
            }
        }

        /// <summary>
        /// Gets or sets the absolute drop tolerance which indicates below what absolute value
        /// an entry will be dropped from the matrix. The standard value is 0.0001.
        /// </summary>
        /// <remarks>
        /// <para>
        /// The values should always be positive and can be larger than 1.0. A low value will
        /// keep more small numbers in the preconditioner matrix. A high value will remove 
        /// more small numbers from the preconditioner matrix.
        /// </para>
        /// <para>
        /// Note that any changes to the <b>DropTolerance</b> after creating the preconditioner
        /// will invalidate the created preconditioner and will require a re-initialization of
        /// the preconditioner.
        /// </para>
        /// </remarks>
        /// <exception cref="ArgumentOutOfRangeException">Thrown if a negative value is provided.</exception>
        public double DropTolerance
        {
            get
            {
                return m_DropTolerance;
            }
            set
            {
                if (value < 0)
                {
                    throw new ArgumentOutOfRangeException("Value");
                }

                m_DropTolerance = value;
            }
        }

        /// <summary>
        /// Gets or sets the pivot tolerance which indicates at what level pivoting will
        /// take place. The standard value is 0.0 which means pivoting will never take place.
        /// </summary>
        /// <remarks>
        /// <para>
        /// The pivot tolerance is used to calculate if pivoting is necessary. Pivoting
        /// will take place if any of the values in a row is bigger than the 
        /// diagonal value of that row divided by the pivot tolerance, i.e. pivoting
        /// will take place if <b>row(i,j) > row(i,i) / PivotTolerance</b> for
        /// any <b>j</b> that is not equal to <b>i</b>.
        /// </para>
        /// <para>
        /// Note that any changes to the <b>PivotTolerance</b> after creating the preconditioner
        /// will invalidate the created preconditioner and will require a re-initialization of
        /// the preconditioner.
        /// </para>
        /// </remarks>
        /// <exception cref="ArgumentOutOfRangeException">Thrown if a negative value is provided.</exception>
        public double PivotTolerance
        {
            get
            {
                return m_PivotTolerance;
            }
            set
            {
                if (value < 0)
                {
                    throw new ArgumentOutOfRangeException("Value");
                }

                m_PivotTolerance = value;
            }
        }
        
        /// <summary>
        /// Returns the upper triagonal matrix that was created during the LU decomposition.
        /// </summary>
        /// <remarks>
        /// This method is used for debugging purposes only and 
        /// should normally not be used.
        /// </remarks>
        /// <returns>A new matrix containing the upper triagonal elements.</returns>
        internal Matrix UpperTriangle()
        {
            return m_Upper.Clone();
        }

        /// <summary>
        /// Returns the lower triagonal matrix that was created during the LU decomposition.
        /// </summary>
        /// <remarks>
        /// This method is used for debugging purposes only and 
        /// should normally not be used.
        /// </remarks>
        /// <returns>A new matrix containing the lower triagonal elements.</returns>
        internal Matrix LowerTriangle()
        {
            return m_Lower.Clone();
        }

        /// <summary>
        /// Returns the pivot array. This array is not needed for normal use because
        /// the preconditioner will return the solution vector values in the proper order.
        /// </summary>
        /// <remarks>
        /// This method is used for debugging purposes only and 
        /// should normally not be used.
        /// </remarks>
        /// <returns>The pivot array.</returns>
        internal int[] Pivots()
        {
            int[] result = new int[m_Pivots.Length];
            for (int i = 0; i < m_Pivots.Length; i++)
            {
                result[i] = m_Pivots[i];
            }

            return result;
        }

        /// <summary>
        /// Initializes the preconditioner and loads the internal data structures.
        /// </summary>
        /// <param name="matrix">
        /// The <see cref="Matrix"/> upon which this preconditioner is based. Note that the 
        /// method takes a general matrix type. However internally the data is stored 
        /// as a sparse matrix. Therefore it is not recommended to pass a dense matrix.
        /// </param>
        /// <exception cref="ArgumentNullException">
        /// If <paramref name="matrix"/> is <b>null</b>.
        /// </exception>
        /// <exception cref="MatrixNotSquareException">
        /// If <paramref name="matrix"/> is not a square matrix.
        /// </exception>
        public void Initialize(Matrix matrix)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("coefficientMatrix");
            }

            if (matrix.Rows != matrix.Columns)
            {
                throw new MatrixNotSquareException();
            }

            SparseMatrix sparseMatrix = (matrix is SparseMatrix) ? matrix as SparseMatrix : new SparseMatrix(matrix);

            // The creation of the preconditioner follows the following algorithm.
            //spaceLeft = lfilNnz * nnz(A)
            //for i = 1, .. , n
            //{
            //    w = a(i,*)
            //    for j = 1, .. , i - 1
            //    {
            //        if (w(j) != 0)
            //        {
            //            w(j) = w(j) / a(j,j)
            //            if (w(j) < dropTol)
            //            {
            //                w(j) = 0;
            //            }
            //            if (w(j) != 0)
            //            {
            //                w = w - w(j) * U(j,*)
            //            }
            //        }
            //    }
            //
            //    for j = i, .. ,n
            //    {
            //        if w(j) <= dropTol * ||A(i,*)||
            //        {
            //            w(j) = 0
            //        }
            //    }
            //
            //    spaceRow = spaceLeft / (n - i + 1) // Determine the space for this row
            //    lfil = spaceRow / 2  // space for this row of L
            //    l(i,j) = w(j) for j = 1, .. , i -1 // only the largest lfil elements
            //
            //    lfil = spaceRow - nnz(L(i,:))  // space for this row of U
            //    u(i,j) = w(j) for j = i, .. , n // only the largest lfil - 1 elements
            //    w = 0
            //    
            //    if max(U(i,i + 1: n)) > U(i,i) / pivTol then // pivot if necessary
            //    {
            //        pivot by swapping the max and the diagonal entries
            //        Update L, U
            //        Update P
            //    }
            //    spaceLeft = spaceLeft - nnz(L(i,:)) - nnz(U(i,:))
            //}

            m_Lower = new SparseMatrix(sparseMatrix.Rows);
            // Create the upper triangular matrix and copy the values
            m_Upper = new SparseMatrix(sparseMatrix.Rows);
            // Create the pivot array
            m_Pivots = new int[sparseMatrix.Rows];
            for (int i = 0; i < m_Pivots.Length; i++)
            {
                m_Pivots[i] = i;
            }
            Vector workVector = new DenseVector(sparseMatrix.Rows);
            Vector rowVector = new DenseVector(sparseMatrix.Columns);
            int[] indexSorting = new int[sparseMatrix.Rows];
            // spaceLeft = lfilNnz * nnz(A)
            int spaceLeft = (int)m_FillLevel * sparseMatrix.NonZeros;
            // for i = 1, .. , n
            for (int i = 0; i < sparseMatrix.Rows; i++)
            {
                // w = a(i,*)
                sparseMatrix.GetRow(i, workVector);
                // pivot the row
                PivotRow(workVector);
                double vectorNorm = workVector.InfinityNorm();
                // for j = 1, .. , i - 1
                for (int j = 0; j < i; j++)
                {
                    // if (w(j) != 0)
                    // {
                    //     w(j) = w(j) / a(j,j)
                    //     if (w(j) < dropTol)
                    //     {
                    //         w(j) = 0;
                    //     }
                    //     if (w(j) != 0)
                    //     {
                    //         w = w - w(j) * U(j,*)
                    //     }
                    if (workVector[j] != 0.0)
                    {
                        // Calculate the multiplication factors that go into 
                        // the L matrix
                        workVector[j] = workVector[j] / m_Upper[j, j];
                        if (System.Math.Abs(workVector[j]) < m_DropTolerance)
                        {
                            workVector[j] = 0.0;
                        }
                        // Calculate the addition factor
                        if (workVector[j] != 0.0)
                        {
                            // vector update all in one go
                            m_Upper.GetRow(j, rowVector);
                            // zero out columnVector[k] because we don't need that
                            // one anymore for k = 0 to k = j
                            for (int k = 0; k <= j; k++)
                            {
                                rowVector[k] = 0.0;
                            }
                            rowVector.Multiply(workVector[j]);
                            workVector.Subtract(rowVector);
                        }
                    }
                }

                //    for j = i, .. ,n
                for (int j = i; j < sparseMatrix.Rows; j++)
                {
                    // if w(j) <= dropTol * ||A(i,*)||
                    // {
                    //     w(j) = 0
                    // }
                    if (System.Math.Abs(workVector[j]) <= m_DropTolerance * vectorNorm)
                    {
                        workVector[j] = 0.0;
                    }
                }

                // spaceRow = spaceLeft / (n - i + 1) // Determine the space for this row
                int spaceRow = spaceLeft / (sparseMatrix.Rows - i + 1);
                // lfil = spaceRow / 2  // space for this row of L
                int fillLevel = spaceRow / 2;
                FindLargestItems(0, i - 1, indexSorting, workVector);
                
                // l(i,j) = w(j) for j = 1, .. , i -1 // only the largest lfil elements
                int lowerNonZeroCount = 0;
                int count = 0;
                for (int j = 0; j < i; j++)
                {
                    if ((count > fillLevel) || (indexSorting[j] == -1))
                    {
                        break;
                    }
                    m_Lower[i, indexSorting[j]] = workVector[indexSorting[j]];
                    count += 1;
                    lowerNonZeroCount += 1;
                }

                FindLargestItems(i + 1, sparseMatrix.Rows - 1, indexSorting, workVector);
                // lfil = spaceRow - nnz(L(i,:))  // space for this row of U
                fillLevel = spaceRow - lowerNonZeroCount;
                // u(i,j) = w(j) for j = i + 1, .. , n // only the largest lfil - 1 elements
                int upperNonZeroCount = 0;
                count = 0;
                for (int j = 0; j < sparseMatrix.Rows - i; j++)
                {
                    if ((count > fillLevel - 1) || (indexSorting[j] == -1))
                    {
                        break;
                    }
                    m_Upper[i, indexSorting[j]] = workVector[indexSorting[j]];
                    count += 1;
                    upperNonZeroCount += 1;
                }

                // Simply copy the diagonal element. Next step is to 
                // see if we pivot
                m_Upper[i, i] = workVector[i];
                // if max(U(i,i + 1: n)) > U(i,i) / pivTol then // pivot if necessary
                // {
                //     pivot by swapping the max and the diagonal entries
                //     Update L, U
                //     Update P
                // }

                // Check if we really need to pivot. If (i+1) >=(mCoefficientMatrix.Rows -1) then
                // we are working on the last row. That means that there is only one number
                // And pivoting is useless. Also the indexSorting array will only contain
                // -1 values.
                if ((i + 1) < (sparseMatrix.Rows - 1))
                {
                    if (System.Math.Abs(workVector[i]) < m_PivotTolerance * System.Math.Abs(workVector[indexSorting[0]]))
                    {
                        // swap columns of u (which holds the values of A in the
                        // sections that haven't been partitioned yet.
                        SwapColumns(m_Upper, i, indexSorting[0]);
                        //  Update P
                        int temp = m_Pivots[i];
                        m_Pivots[i] = m_Pivots[indexSorting[0]];
                        m_Pivots[indexSorting[0]] = temp;
                    }
                }
                // spaceLeft = spaceLeft - nnz(L(i,:)) - nnz(U(i,:))
                spaceLeft -= (lowerNonZeroCount + upperNonZeroCount);
            }

            for (int i = 0; i < m_Lower.Rows; i++)
            {
                m_Lower[i, i] = 1.0;
            }
        }

        private void PivotRow(Vector row)
        {
            Dictionary<int, int> knownPivots = new Dictionary<int, int>();
            // pivot the row
            for (int i = 0; i < row.Count; i++)
            {
                if ((m_Pivots[i] != i) && (!PivotMapFound(knownPivots, i)))
                {
                    // store the pivots in the hashtable
                    knownPivots.Add(m_Pivots[i], i);

                    double t = row[i];
                    row[i] = row[m_Pivots[i]];
                    row[m_Pivots[i]] = t;
                }
            }
        }

        private bool PivotMapFound(Dictionary<int, int> knownPivots, int currentItem)
        {
            if (knownPivots.ContainsKey(m_Pivots[currentItem]))
            {
                if (knownPivots[m_Pivots[currentItem]].Equals(currentItem))
                {
                    return true;
                }
            }

            if ((knownPivots.ContainsKey(currentItem)))
            {
                if (knownPivots[currentItem].Equals(m_Pivots[currentItem]))
                {
                    return true;
                }
            }

            return false;
        }

        private static void SwapColumns(Matrix matrix, int firstColumn, int secondColumn)
        {
            for (int i = 0; i < matrix.Rows; i++)
            {
                double temp = matrix[i, firstColumn];
                matrix[i, firstColumn] = matrix[i, secondColumn];
                matrix[i, secondColumn] = temp;
            }
        }

        private static void FindLargestItems(int lowerBound, int upperBound, int[] sortedIndices, Vector values)
        {
            // Copy the indices for the values into the array
            for (int i = 0; i < upperBound + 1 - lowerBound; i++)
            {
                sortedIndices[i] = lowerBound + i;
            }
            for (int i = upperBound + 1 - lowerBound; i < sortedIndices.Length; i++)
            {
                sortedIndices[i] = -1;
            }
            // Sort the first set of items.
            // Sorting starts at index 0 because the index array
            // starts at zero
            // and ends at index upperBound - lowerBound
            IlutpElementSorter.SortDoubleIndicesDecreasing(0, upperBound - lowerBound, sortedIndices, values);
        }

        /// <summary>
        /// Approximates the solution to the matrix equation <b>Ax = b</b>.
        /// </summary>
        /// <param name="rhs">The right hand side vector.</param>
        /// <returns>The left hand side vector.</returns>
        public Vector Approximate(Vector rhs)
        {
            if (rhs == null)
            {
                throw new ArgumentNullException("rhs");
            }

            if (m_Upper == null)
            {
                throw new MissingMatrixException();
            }

            if (rhs.Count != m_Upper.Columns)
            {
                throw new NotConformableException();
            }

            Vector result = new DenseVector(rhs.Count);
            Approximate(rhs, result);
            return result;
        }

        /// <summary>
        /// Approximates the solution to the matrix equation <b>Ax = b</b>.
        /// </summary>
        /// <param name="rhs">The right hand side vector.</param>
        /// <param name="lhs">The left hand side vector. Also known as the result vector.</param>
        public void Approximate(Vector rhs, Vector lhs)
        {
            if (rhs == null)
            {
                throw new ArgumentNullException("rhs");
            }

            if (lhs == null)
            {
                throw new ArgumentNullException("lhs");
            }

            if (m_Upper == null)
            {
                throw new MissingMatrixException();
            }

            if ((lhs.Count != rhs.Count) || (lhs.Count != m_Upper.Rows))
            {
                throw new NotConformableException();
            }
            
            // Solve equation here
            // Pivot(vector, result);
            // Solve L*Y = B(piv,:)
            Vector rowValues = new DenseVector(m_Lower.Rows);
            for (int i = 0; i < m_Lower.Rows; i++)
            {
                m_Lower.GetRow(i, rowValues);

                double sum = 0.0;
                for (int j = 0; j < i; j++)
                {
                    sum += rowValues[j] * lhs[j];
                }

                lhs[i] = (rhs[i] - sum);
            }

            // Solve U*X = Y;
            for (int i = m_Upper.Rows - 1; i > -1; i--)
            {
                m_Upper.GetRow(i, rowValues);

                double sum = 0.0;
                for (int j = m_Upper.Rows - 1; j > i; j--)
                {
                    sum += rowValues[j] * lhs[j];
                }

                lhs[i] = 1 / rowValues[i] * (lhs[i] - sum);
            }

            // We have a column pivot so we only need to pivot the
            // end result not the incoming right hand side vector
            Vector temp = lhs.Clone();
            Pivot(temp, lhs);
        }

        private void Pivot(Vector vector, Vector result)
        {
            for (int i = 0; i < m_Pivots.Length; i++)
            {
                result[i] = vector[m_Pivots[i]];
            }
        }
    }
}
